1 Study description

Microbial and faunal assemblages associated with the tubeworm Ridgeia piscesae were collected from thirteen tubeworm bushes in the Endeavour Hydrothermal Vent Marine Protected Area and Middle Valley aboard the CCGS Tully using the remotely-operated vehicle ROPOS in August of 2016. Macro- and meiofaunal species were identified and enumerated from 9 assemblage (“grab”) samples. Microbial community DNA was extracted from biofilm/detrital matter associated with all 13 assemblages. Microbial biomass was separated into two size fractions (20-64µm, 0.2-20µm) by serial filtration and two method replicates for each size fraction comparing biomass removal methods were included from assemblages EMw1 and EMw6. Microbial DNA was extracted from 7 diffuse fluid and 3 background fluid samples for a total of 40 DNA extractions.

2 Data Import

2.1 Faunal species data

The faunal dataset contains counts of individuals for 41 macrofaunal taxa (1mm sieve) and 25 mieofaunal taxa (64µm sieve; inlcudes juveniles of 11 macrofauna). Large samples were split, counted and scaled up to full sample numbers. Details of counted splits are in Supplementary Tables 1 and 2. Data are provided as counts for meiofauna and relative abundance for macrofauna due to planned inclusion of macrofaunal species counts in another manuscript. Please contact R. Boschen-Rose () for macrofaunal species counts.

# Macrofauna species counts (available by request)
macro9 <- read.csv("Macro9counts.csv", header = TRUE, row.names = 1)
macro9 <- as.data.frame(t(macro9))

# Meiofauna species counts
meio9 <- read.csv("Meio9counts.csv", header = TRUE, row.names = 1)
meio9 <- as.data.frame(t(meio9))

2.2 Microbial OTU data

The microbial dataset contains OTU counts from reads produced on Illumina MiSeq runs using 16S/18S rRNA gene primers for Bacteria, Archaea, and microeukaryotes. Size fractions and method replicates are combined here. See “SamplingBiasTests” for justification.

# Bacterial OTU counts
codaBact23nonsingleton <- read.csv("BacterialOTUs.csv", header = TRUE, 
    row.names = 1)
codaBact23nonsingleton <- as.data.frame(t(codaBact23nonsingleton))

# Eukaryal OTU counts
codaEuk23nonsingleton <- read.csv("EukaryalOTUs.csv", header = TRUE, row.names = 1)
codaEuk23nonsingleton <- as.data.frame(t(codaEuk23nonsingleton))

# Archaeal OTU counts (ECw11 not included- too few reads)
codaArch22nonsingleton <- read.csv("ArchaealOTUs.csv", header = TRUE, row.names = 1)
codaArch22nonsingleton <- as.data.frame(t(codaArch22nonsingleton))

Files with size fractions and replicates kept separate (40 individual DNA extracts). Singleton OTUs have been removed.

# Bacteria
Bact40 <- read.csv("Bact40nonsingleton.csv", header = TRUE, row.names = 1)

# Eukarya
Euk40 <- read.csv("Euk40nonsingleton.csv", header = TRUE, row.names = 1)

# Archaea (two extracts of ECw11 removed- too few reads)
Arch38 <- read.csv("Arch38nonsingleton.csv", header = TRUE, row.names = 1)

QPCR-balanced microbial community. Three microbial domains were balanced to whole community relative abundances based on quantitative PCR results and multiplied by 105 to scale to whole number pseudo counts.

balMic <- read.csv("QPCR_balanced_microbes.csv", header = TRUE, row.names = 1)
balMic <- as.data.frame(t(balMic))

2.3 Metadata

Sample_type codes: f=fluids, B=background, bk=cell removal through agitation in bucket of artificial seawater, d=direct removal of cells in biofilms. T_category for diffuse fluids is determined by the basal temperature of the associated assemblage, NOT the diffuse fluid temperature itself.

metadata23 <- read.csv("metadata23.csv", header = TRUE, row.names = 1)
metadata23
##           Source Sample_type     Vent_field Temperature T_category Substratum
## EMw1   Tubeworms          bk Main_Endeavour        33.2       high    sulfide
## EMw2   Tubeworms           d Main_Endeavour         5.3        low    sulfide
## EMw3   Tubeworms          bk Main_Endeavour         3.4        low     basalt
## EMw4   Tubeworms           d Main_Endeavour        37.1       high    sulfide
## EMw5   Tubeworms          bk Main_Endeavour        13.4        low    sulfide
## EMw6   Tubeworms          bk Main_Endeavour         2.8        low     basalt
## EMw7   Tubeworms           d Main_Endeavour        33.0       high    sulfide
## EMw8   Tubeworms           d Main_Endeavour        37.4       high    sulfide
## ECw9   Tubeworms          bk       Clam_Bed        15.0        low    sulfide
## ECw10  Tubeworms          bk       Clam_Bed        31.0       high    sulfide
## ECw11  Tubeworms           d       Clam_Bed         5.2        low     basalt
## MVw12  Tubeworms          bk  Middle_Valley        27.4       high    sulfide
## MVw13  Tubeworms          bk  Middle_Valley        21.0        low    sulfide
## EMf1     Diffuse           f Main_Endeavour         5.0        low    sulfide
## EMf2     Diffuse           f Main_Endeavour         4.6        low    sulfide
## EMf3     Diffuse           f Main_Endeavour         4.1        low     basalt
## EMf7     Diffuse           f Main_Endeavour        12.0        low    sulfide
## ECf9     Diffuse           f       Clam_Bed         2.6        low    sulfide
## ECf11    Diffuse           f       Clam_Bed         6.1        low     basalt
## MVf12    Diffuse           f  Middle_Valley        23.0        low    sulfide
## bkEC  Background           B       Clam_Bed         2.0        low       none
## bkEM  Background           B Main_Endeavour         2.0        low       none
## bkMV  Background           B  Middle_Valley         2.0        low       none

2.4 Load necessary R packages

library(vegan)
library(reshape2)
library(ggplot2)
library(ggrepel)
library(ggpubr)
library(funrar)
library(compositions)
library(zCompositions)
library(propr)
library(ALDEx2)


3 Sample diversity

Diversity was calculated using the Inverse Simpson metric. For microbial domains from tubeworm-associated samples, diversity is calculated individually on all extracts and then with size fractions and method replicates combined. Wilcoxon Rank Sum tests are applied to test for significant differences in diversity levels.

3.1 Macrofauna

MacroDiv <- diversity(macro9, index = "invsimpson", MARGIN = 1)
MacroDiv
##     EMw3    ECw11     EMw1     EMw4     EMw8    MVw12     EMw5    ECw10 
## 2.357427 2.626470 3.031563 3.046310 3.043524 2.297812 1.266767 1.940417 
##    MVw13 
## 9.226349

Wilcox test between highT and lowT diversity levels.

MAChighTdiv <- c(3.03, 3.04, 3.04, 2.3, 1.94)
MAClowTdiv <- c(2.36, 2.62, 1.27, 9.23)

wilcox.test(MAClowTdiv, MAChighTdiv)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  MAClowTdiv and MAChighTdiv
## W = 9, p-value = 0.9021
## alternative hypothesis: true location shift is not equal to 0

Macrofaunal diversity IS NOT significantly different between highT and lowT grabs.

3.2 Meiofauna

MeioDiv <- diversity(meio9, index = "invsimpson", MARGIN = 1)
MeioDiv
##     EMw3    ECw11     EMw1     EMw4     EMw8    MVw12     EMw5    ECw10 
## 4.784573 4.867552 1.242057 1.078020 1.012415 1.000000 3.020960 1.380738 
##    MVw13 
## 1.543465

Wilcox test between highT and lowT diversity levels.

MEhighTdiv <- c(1.24, 1.07, 1.01, 1, 1.38)
MElowTdiv <- c(4.78, 4.86, 3.02, 1.54)

wilcox.test(MEhighTdiv, MElowTdiv)
## 
##  Wilcoxon rank sum exact test
## 
## data:  MEhighTdiv and MElowTdiv
## W = 0, p-value = 0.01587
## alternative hypothesis: true location shift is not equal to 0

Meiofaunal diversity IS significantly different between highT and lowT grabs.

3.3 Microbes (40 extracts)

For individual DNA extracts, sample naming conventions for tubeworm associated extracts are as follows:
Sample name followed by microbial community removal method (bk= agitation in bucket of artificial seawater, d= direct removal of attached biofilm) and size fraction (m= micro 20-64µm, pn= pico/nano 0.2-20µm).
EXAMPLE: EMw3bkm –> sample EMw3, bucket agitation method, micro size fraction.

3.3.1 Bacteria

Bact40nonsingDiv <- diversity(Bact40, index = "invsimpson", MARGIN = 1)
Bact40nonsingDiv
##    EMw3bkm   EMw3bkpn       EMf3    EMw6bkm   EMw6bkpn     EMw6dm    EMw6dpn 
##  13.877373 132.019034   3.356216 190.791269 129.672086  96.763474 180.551378 
##    ECw11dm   ECw11dpn      ECf11    EMw1bkm   EMw1bkpn     EMw1dm    EMw1dpn 
##  18.147660  12.664641  13.794946  69.317711  67.080719  13.660772  22.860723 
##       EMf1     EMw4dm    EMw4dpn     EMw7dm    EMw7dpn       EMf7     EMw8dm 
##   6.725429   5.003208   5.384283  30.205731  16.667821   8.739932  26.253335 
##    EMw8dpn   MVw12bkm  MVw12bkpn      MVf12     EMw2dm    EMw2dpn       EMf2 
##  30.179900  11.177898  13.907539   6.783713 114.013041  61.603672   6.001158 
##    EMw5bkm   EMw5bkpn    ECw9bkm   ECw9bkpn       ECf9   ECw10bkm  ECw10bkpn 
## 108.581991  13.380956  38.195485 125.157585  10.439104  10.566329  15.878173 
##   MVw13bkm  MVw13bkpn       bkCB       bkEM       bkMV 
##  30.641508  39.521918   9.444129  13.949291  16.312045
BhighTdiv <- c(11.18, 10.57, 69.32, 30.21, 13.66, 5, 26.25, 13.91, 15.88, 
    67.08, 16.67, 22.86, 5.38, 30.18)
BlowTdiv <- c(190.79, 13.88, 108.58, 38.2, 30.64, 96.76, 18.15, 114.01, 
    129.67, 132.02, 13.38, 125.16, 39.52, 180.55, 12.66, 61.6)
BfluidsDiv <- c(9.44, 13.95, 16.31, 3.36, 13.79, 6.73, 8.74, 6.78, 6, 10.44)

wilcox.test(BhighTdiv, BlowTdiv)
## 
##  Wilcoxon rank sum exact test
## 
## data:  BhighTdiv and BlowTdiv
## W = 43, p-value = 0.003321
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(BfluidsDiv, BhighTdiv)
## 
##  Wilcoxon rank sum exact test
## 
## data:  BfluidsDiv and BhighTdiv
## W = 30, p-value = 0.01855
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(BfluidsDiv, BlowTdiv)
## 
##  Wilcoxon rank sum exact test
## 
## data:  BfluidsDiv and BlowTdiv
## W = 8, p-value = 2.523e-05
## alternative hypothesis: true location shift is not equal to 0

Bacterial diversity IS signficantly different between highT and lowT grabs and between fluids and both highT and lowT grabs.

3.3.2 Eukarya

Euk40nonsingDiv <- diversity(Euk40, index = "invsimpson", MARGIN = 1)
Euk40nonsingDiv
##   EMw3bkm  EMw3bkpn      EMf3   EMw6bkm  EMw6bkpn    EMw6dm   EMw6dpn   ECw11dm 
## 24.935378  7.484512 20.319558 11.688583  2.682649  2.561250  3.220499  9.205386 
##  ECw11dpn     ECf11   EMw1bkm  EMw1bkpn    EMw1dm   EMw1dpn      EMf1    EMw4dm 
##  3.037183  7.435008  8.593577  2.427560  1.325924  1.225294  5.474317  4.348228 
##   EMw4dpn    EMw7dm   EMw7dpn      EMf7    EMw8dm   EMw8dpn MVw12bkpn   MVw12bm 
##  3.277309  2.936878  3.693987 14.214942  2.701855  5.718646 11.698233  5.852115 
##     MVf12    EMw2dm   EMw2dpn      EMf2   EMw5bkm  EMw5bkpn   ECw9bkm  ECw9bkpn 
##  8.343386  8.200491  8.587284 16.754450 27.093450 40.638687  8.334630 17.235464 
##      ECf9  ECw10bkm ECw10bkpn  MVw13bkm MVw13bkpn      bkCB      bkEM      bkMV 
##  1.751818  4.597617 13.682620  6.739444  8.104348  2.013199  1.829072 33.575884
EhighTdiv <- c(11.7, 4.6, 2.94, 8.59, 1.33, 4.35, 2.7, 5.85, 13.68, 3.69, 
    2.43, 1.23, 3.28, 5.72)
ElowTdiv <- c(11.69, 2.56, 24.94, 9.21, 8.2, 27.09, 8.33, 6.74, 2.68, 3.22, 
    7.48, 3.04, 8.59, 40.64, 17.24, 8.1)
EfluidsDiv <- c(2.01, 1.83, 33.58, 20.32, 7.44, 5.47, 14.21, 8.34, 16.75, 
    1.75)

wilcox.test(EhighTdiv, ElowTdiv)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  EhighTdiv and ElowTdiv
## W = 61.5, p-value = 0.03764
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(EfluidsDiv, EhighTdiv)
## 
##  Wilcoxon rank sum exact test
## 
## data:  EfluidsDiv and EhighTdiv
## W = 93, p-value = 0.1915
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(ElowTdiv, EfluidsDiv)
## 
##  Wilcoxon rank sum exact test
## 
## data:  ElowTdiv and EfluidsDiv
## W = 90, p-value = 0.6227
## alternative hypothesis: true location shift is not equal to 0

Eukaryal diversity IS signficantly different between highT and lowT grabs, but not between fluids and either highT or lowT grabs.

3.3.3 Archaea

Arch40nonsingDiv <- diversity(Arch38, index = "invsimpson", MARGIN = 1)
Arch40nonsingDiv
##   EMw3bkm  EMw3bkpn    EMf3pn   EMw6bkm  EMw6bkpn    EMw6dm   EMw6dpn   ECf11pn 
##  1.020409  1.034489  4.710203  1.042860  1.032802  1.182467  1.012818  3.541635 
##   EMw1bkm  EMw1bkpn    EMw1dm   EMw1dpn    EMf1pn    EMw4dm   EMw4dpn    EMw7dm 
##  1.121875  1.619707  6.671790  6.443250  3.614836  2.096764  6.405379  6.050011 
##   EMw7dpn    EMf7pn    EMw8dm   EMw8dpn  MVw12bkm MVw12bkpn   MVf12pn    EMw2dm 
##  7.688165  3.067253  5.994080  4.717821  3.940932  7.645768  4.493723  2.341680 
##   EMw2dpn    EMf2pn   EMw5bkm  EMw5bkpn   ECw9bkm  ECw9bkpn    ECf9pn  ECw10bkm 
##  1.030438  4.646234  1.018484  1.014421  4.146493  2.351466  4.397208  8.000000 
## ECw10bkpn  MVw13bkm MVw13bkpn    bkECpn    bkEMpn    bkMVpn 
##  9.403731  2.860098  2.608518  3.242648  2.936804  5.076730
AhighTdiv <- c(3.94, 8, 6.05, 1.12, 6.67, 2.1, 5.99, 7.65, 9.4, 7.69, 1.62, 
    6.44, 6.41, 4.72)
AlowTdiv <- c(1.04, 1.18, 1.02, 2.34, 1.02, 4.15, 2.86, 1.03, 1.01, 1.03, 
    1.03, 1.01, 2.35, 2.61)
AfluidsDiv <- c(3.24, 2.94, 5.08, 4.71, 3.54, 3.61, 3.07, 4.49, 4.65, 4.4)

wilcox.test(AlowTdiv, AhighTdiv)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  AlowTdiv and AhighTdiv
## W = 17, p-value = 0.0002141
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(AfluidsDiv, AhighTdiv)
## 
##  Wilcoxon rank sum exact test
## 
## data:  AfluidsDiv and AhighTdiv
## W = 36, p-value = 0.0484
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(AlowTdiv, AfluidsDiv)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  AlowTdiv and AfluidsDiv
## W = 5, p-value = 0.0001558
## alternative hypothesis: true location shift is not equal to 0

Archaeal diversity IS signficantly different between highT and lowT grabs and between fluids and both highT and lowT grabs.

3.4 Diversity boxplots

(Plotting for all 40 extracts)

# Combine diversity values into one dataframe
InvSimpson <- c(BhighTdiv, BlowTdiv, BfluidsDiv, EhighTdiv, ElowTdiv, EfluidsDiv, 
    AhighTdiv, AlowTdiv, AfluidsDiv, MAChighTdiv, MAClowTdiv, MEhighTdiv, 
    MElowTdiv)
domain <- c(rep("Bacteria", 40), rep("Eukarya", 40), rep("Archaea", 38), 
    rep("Macrofauna", 9), rep("Meiofauna", 9))
Sample_type <- c(rep(c(rep("highT", 14), rep("lowT", 16), rep("fluids", 
    10)), 2), rep("highT", 14), rep("lowT", 14), rep("fluids", 10), rep(c(rep("highT", 
    5), rep("lowT", 4)), 2))

Div4Plots <- data.frame(InvSimpson, Sample_type, domain)

# Plot them
ggplot(Div4Plots, aes(domain, InvSimpson, fill = Sample_type)) + geom_boxplot(varwidth = FALSE) + 
    scale_fill_discrete(breaks = c("fluids", "highT", "lowT"), labels = c("Fluids", 
        "HighT", "LowT"), name = "Sample Type") + labs(x = NULL, y = "Diversity (Inverse Simpson)") + 
    facet_wrap(. ~ domain, scales = "free") + theme(axis.text.x = element_blank(), 
    axis.title.x = element_blank(), axis.text.y = element_text(size = 12), 
    axis.title.y = element_text(size = 14))


3.5 Microbes (23 samples)

For these calculations, size fractions and method replicates are combined.

3.5.1 Bacteria

NOTE: These diversity calculations are based on subsampled datasets (using the smallest read count by domain).

rowSums(codaBact23nonsingleton)
##   EMf3  ECf11   EMf1   EMf7  MVf12   EMf2   ECf9   bkEC   bkEM   bkMV   EMw3 
##  34535  20711  28350  30036  18755  34250  24775  21278  34180  26397  44531 
##   EMw6  ECw11   EMw1   EMw4   EMw7   EMw8  MVw12   EMw2   EMw5   ECw9  ECw10 
##  83006  32677 112883  58955  60625  64972  45872  66446  30031  49338  51516 
##  MVw13 
##  24247
subBact23 <- as.data.frame(rrarefy(codaBact23nonsingleton, sample = 18755))
subBact23Div <- diversity(subBact23, index = "invsimpson", MARGIN = 1)
subBact23Div
##       EMf3      ECf11       EMf1       EMf7      MVf12       EMf2       ECf9 
##   3.380745  13.753857   6.662416   8.823526   6.783713   5.940150  10.509222 
##       bkEC       bkEM       bkMV       EMw3       EMw6      ECw11       EMw1 
##   9.462509  14.124186  16.546602  36.344101 217.909681  16.925666  47.139049 
##       EMw4       EMw7       EMw8      MVw12       EMw2       EMw5       ECw9 
##   5.394489  25.912545  27.535837  13.170838  94.589946  95.946261  73.961936 
##      ECw10      MVw13 
##  12.918118  37.545077

3.5.2 Eukarya

subEuk23 <- as.data.frame(rrarefy(codaEuk23nonsingleton, sample = 13172))
subEuk23Div <- diversity(subEuk23, index = "invsimpson", MARGIN = 1)
subEuk23Div
##      EMf3     ECf11      EMf1      EMf7     MVf12      EMf2      ECf9      bkEC 
## 21.132670  7.371347  5.328299 14.430674  8.329511 17.015357  1.772109  1.984034 
##      bkEM      bkMV      EMw3      EMw6     ECw11      EMw1      EMw4      EMw7 
##  1.863651 34.021920  8.507661  5.353393  8.783340  4.391468  6.282467  4.685926 
##      EMw8     MVw12      EMw2      EMw5      ECw9     ECw10     MVw13 
##  6.034257  8.896173 10.016093 42.089294 15.393216 10.769617  7.604550

3.5.3 Archaea

subArch22 <- as.data.frame(rrarefy(codaArch22nonsingleton, sample = 25064))
subArch22Div <- diversity(subArch22, index = "invsimpson", MARGIN = 1)
subArch22Div
##     EMf3    ECf11     EMf1     EMf7    MVf12     EMf2     ECf9     bkEC 
## 4.711096 3.512507 3.605891 3.024852 4.461340 4.650887 4.486183 3.278287 
##     bkEM     bkMV     EMw3     EMw6     EMw1     EMw4     EMw7     EMw8 
## 2.907239 5.035000 1.028663 1.025469 1.805728 5.709097 6.991795 4.897325 
##    MVw12     EMw2     EMw5     ECw9    ECw10    MVw13 
## 5.407972 1.031313 1.014593 2.875343 9.411689 2.731627


4 Comparing sample compositions

4.1 Compositional data approach

Data from high-throughput sequencing platforms, which have an arbitrary sum imposed by the instrument, contain only information on ratios between taxa and are therefore considered compositional in nature. Data analysis followed the compositional approach recommended by Gloor and Reid (2016) and Quinn et al (2019). Zero counts were replaced by an imputed value using the count zero multiplicative method in the zCompositions R package (Martín-Fernández et al 2011, Palarea-Albaladejo and Martin-Fernandez 2015), and data were then transformed by centered log-ratio (Aitchison 1986). Faunal data followed the same approach because full assemblages were not collected.

4.2 Macrofauna

Nonmetric Multidimensional Scaling to visualize the relationships between samples based on species abundances.

# Package 'zCompositions' for imputation of zeros. Output='p-counts'
# returns pseudo-counts (preserves ratios in the original data after
# conversion of zeros to non-zeros).
mac9_cmultRepl <- cmultRepl(macro9, method = "CZM", output = "p-counts")

# Centered log-ratio transformation
clr_mac9 <- apply(mac9_cmultRepl, 2, function(x) {
    log(x) - mean(log(x))
})

# Package 'compositions' to calculate Aitchison distance matrix metaMDS
# for Nonmetric Multidimensional Scaling
clr_mac9.dist <- dist(clr_mac9)
clr_mac9.NMDS <- metaMDS(clr_mac9.dist, k = 2, try = 100)

# Trim the metadata file down to only the samples with all three size
# classes characterized and re-order the sample rows to match the
# macro9 file.
metadata9 <- metadata23[c(1, 3, 4, 5, 8, 10, 11, 12, 13), ]
metadata9 <- metadata9[match(rownames(macro9), rownames(metadata9)), ]

# Make a dataframe with NMDS points and associated metadata
clr_mac9.NMDS.df <- data.frame(x = clr_mac9.NMDS$points[, 1], y = clr_mac9.NMDS$points[, 
    2], T_category = as.factor(metadata9[, 5]), Vent_field = as.factor(metadata9[, 
    3]), Temperature = metadata9[, 4], Substratum = as.factor(metadata9[, 
    6]))

# Plot it
ggplot(data = clr_mac9.NMDS.df, aes(x = x, y = y, color = Temperature, 
    shape = Substratum)) + geom_point(size = 4) + geom_text_repel(aes(label = rownames(metadata9)), 
    hjust = 0.5, vjust = 1.5, size = 4) + ggtitle("Macrofauna NMDS") + 
    scale_color_gradient(low = "blue", high = "red")


Envfit to see which variables (Vent_field + Substratum + Temperature + T_category) explain the NMDS plot.
Temperature is a signficant predictor, followed by substratum.

clr_mac9_NMDSfit <- envfit(clr_mac9.NMDS ~ ., metadata9)
clr_mac9_NMDSfit
## 
## ***VECTORS
## 
##                NMDS1    NMDS2    r2 Pr(>r)   
## Temperature  0.99151 -0.13002 0.927  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                            NMDS1   NMDS2
## SourceTubeworms           0.0000  0.0000
## Sample_typebk            -0.7698 -0.0033
## Sample_typed              1.5396  0.0066
## Vent_fieldClam_Bed       -0.7717 -0.0030
## Vent_fieldMain_Endeavour  0.6161  0.0021
## Vent_fieldMiddle_Valley  -0.7686 -0.0021
## T_categoryhigh            4.7975 -2.3862
## T_categorylow            -5.9969  2.9827
## Substratumbasalt         -7.7221 -0.0325
## Substratumsulfide         2.2063  0.0093
## 
## Goodness of fit:
##                 r2 Pr(>r)  
## Source      0.0000  1.000  
## Sample_type 0.0171  1.000  
## Vent_field  0.0069  0.988  
## T_category  0.5192  0.012 *
## Substratum  0.2465  0.026 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999


Hierarchical cluster analysis.

clr_mac9.CLUST <- hclust(clr_mac9.dist, method = "average", member = NULL)
plot(clr_mac9.CLUST, hang = -1, main = "Macrofauna")


ANOSIM tests on the clr-transformed data to assess differences based on highT vs lowT, substratum and vent field.

mac9ANOSIMtemp <- with(metadata9, anosim(clr_mac9.dist, T_category))
summary(mac9ANOSIMtemp)
## 
## Call:
## anosim(x = clr_mac9.dist, grouping = T_category) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.7812 
##       Significance: 0.012 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.269 0.325 0.381 0.781 
## 
## Dissimilarity ranks between and within classes:
##         0%   25%  50%   75% 100%  N
## Between 13 17.75 25.5 31.25   36 20
## high     1  3.25  6.0 10.75   23 10
## low      6  8.25 14.5 21.50   28  6
mac9ANOSIMsubstr <- with(metadata9, anosim(clr_mac9.dist, Substratum))
summary(mac9ANOSIMsubstr)
## 
## Call:
## anosim(x = clr_mac9.dist, grouping = Substratum) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.3961 
##       Significance: 0.107 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.396 0.474 0.656 0.656 
## 
## Dissimilarity ranks between and within classes:
##         0%   25%  50%   75% 100%  N
## Between  6 17.25 24.5 29.75   35 14
## basalt   8  8.00  8.0  8.00    8  1
## sulfide  1  7.00 14.0 23.00   36 21
mac9ANOSIMfield <- with(metadata9, anosim(clr_mac9.dist, Vent_field))
summary(mac9ANOSIMfield)
## 
## Call:
## anosim(x = clr_mac9.dist, grouping = Vent_field) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.2847 
##       Significance: 0.091 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.278 0.368 0.458 0.646 
## 
## Dissimilarity ranks between and within classes:
##                0%   25%  50%   75% 100%  N
## Between         4 10.75 21.5 28.25   35 24
## Clam_Bed       32 32.00 32.0 32.00   32  1
## Main_Endeavour  1  3.75 13.5 15.75   25 10
## Middle_Valley  36 36.00 36.0 36.00   36  1

RESULT: Temperature category is a significant factor in macrofaunal diversity variation, substratum and vent field are not.


4.3 Meiofauna


Envfit. Temperature is the only signficant predictor.

## 
## ***VECTORS
## 
##               NMDS1   NMDS2     r2 Pr(>r)   
## Temperature 0.99407 0.10876 0.9072  0.003 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                             NMDS1    NMDS2
## SourceTubeworms            0.0000   0.0000
## Sample_typebk             -1.0118   1.4433
## Sample_typed               2.0236  -2.8866
## Vent_fieldClam_Bed        -3.1965  -3.2965
## Vent_fieldMain_Endeavour  -0.1757  -3.0498
## Vent_fieldMiddle_Valley    3.6357  10.9209
## T_categoryhigh            12.5709  -2.0716
## T_categorylow            -15.7136   2.5895
## Substratumbasalt         -19.1619  -4.5150
## Substratumsulfide          5.4748   1.2900
## 
## Goodness of fit:
##                 r2 Pr(>r)   
## Source      0.0000  1.000   
## Sample_type 0.0218  0.700   
## Vent_field  0.1376  0.738   
## T_category  0.7105  0.008 **
## Substratum  0.3877  0.054 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999


ANOSIM tests on the clr-transformed data (highT vs lowT, substratum and vent field).

## 
## Call:
## anosim(x = clr_meio9.dist, grouping = T_category) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.9562 
##       Significance: 0.004 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.344 0.425 0.444 0.669 
## 
## Dissimilarity ranks between and within classes:
##         0%   25%  50%   75% 100%  N
## Between 15 21.50 26.5 31.25   36 20
## high     1  3.25  5.5 10.00   13 10
## low      8  9.25 12.0 16.25   21  6
## 
## Call:
## anosim(x = clr_meio9.dist, grouping = Substratum) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.5325 
##       Significance: 0.06 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.175 0.532 0.695 0.695 
## 
## Dissimilarity ranks between and within classes:
##         0%  25%  50%  75% 100%  N
## Between  8 19.5 25.5 31.5   36 14
## basalt  10 10.0 10.0 10.0   10  1
## sulfide  1  6.0 15.0 22.0   33 21
## 
## Call:
## anosim(x = clr_meio9.dist, grouping = Vent_field) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: -0.1806 
##       Significance: 0.862 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.285 0.396 0.467 0.507 
## 
## Dissimilarity ranks between and within classes:
##                0%   25%  50%   75% 100%  N
## Between         1  9.75 17.5 23.25   36 24
## Clam_Bed       30 30.00 30.0 30.00   30  1
## Main_Endeavour  3  8.75 25.5 28.75   34 10
## Middle_Valley  15 15.00 15.0 15.00   15  1

RESULT: Temperature category is a significant factor in meiofaunal diversity variation, substratum and vent field are not.


4.4 Microbes (9 fully characterized samples)

4.4.1 Bacteria

codaBact9 <- codaBact23nonsingleton[c(11, 13, 14, 15, 17, 18, 20, 22, 23), 
    ]

# remove any OTU with only 1 read. Drops from 40,389 to 18,570 OTUs.
codaBact9nonsingleton <- codaBact9[, which(colSums(codaBact9) > 1)]


Envfit. Temperature is the only signficant predictor.

## 
## ***VECTORS
## 
##               NMDS1   NMDS2     r2 Pr(>r)   
## Temperature 0.67513 0.73770 0.8323  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                             NMDS1    NMDS2
## SourceTubeworms            0.0000   0.0000
## Sample_typebk             -9.1711 -10.3184
## Sample_typed              18.3422  20.6368
## Vent_fieldClam_Bed        -3.0742  -3.6403
## Vent_fieldMain_Endeavour  16.4566 -17.8509
## Vent_fieldMiddle_Valley  -38.0674  48.2674
## T_categoryhigh            38.6052  21.1261
## T_categorylow            -48.2565 -26.4076
## Substratumbasalt         -46.6324 -58.8916
## Substratumsulfide         13.3235  16.8262
## 
## Goodness of fit:
##                 r2 Pr(>r)  
## Source      0.0000  1.000  
## Sample_type 0.0569  0.734  
## Vent_field  0.1749  0.680  
## T_category  0.3612  0.017 *
## Substratum  0.2405  0.169  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999


ANOSIM tests on the clr-transformed data (Temp category and substratum).

## 
## Call:
## anosim(x = clr_bact9.dist, grouping = Substratum) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.1494 
##       Significance: 0.342 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.448 0.513 0.838 0.838 
## 
## Dissimilarity ranks between and within classes:
##         0%   25%  50%   75% 100%  N
## Between  4 11.75 21.5 29.75   34 14
## basalt  22 22.00 22.0 22.00   22  1
## sulfide  1  9.00 17.0 25.00   36 21
## 
## Call:
## anosim(x = clr_bact9.dist, grouping = T_category) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.325 
##       Significance: 0.03 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.175 0.275 0.325 0.381 
## 
## Dissimilarity ranks between and within classes:
##         0%   25%  50%   75% 100%  N
## Between  5 14.00 19.5 30.25   35 20
## high     1  5.25 17.0 24.75   36 10
## low      4  7.50 11.5 20.00   23  6

RESULT: Temperature category is a weak factor in bacterial diversity variation, substratum is not.


4.4.2 Eukarya


Envfit. No signficant predictors.

## 
## ***VECTORS
## 
##                NMDS1    NMDS2     r2 Pr(>r)  
## Temperature -0.05364  0.99856 0.5827  0.093 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                             NMDS1    NMDS2
## SourceTubeworms            0.0000   0.0000
## Sample_typebk              2.2839  -7.0991
## Sample_typed              -4.5677  14.1982
## Vent_fieldClam_Bed         9.6160   9.8077
## Vent_fieldMain_Endeavour   9.3395 -13.1135
## Vent_fieldMiddle_Valley  -32.9648  22.9760
## T_categoryhigh             1.4697  21.6364
## T_categorylow             -1.8372 -27.0455
## Substratumbasalt          19.4432 -33.6297
## Substratumsulfide         -5.5552   9.6085
## 
## Goodness of fit:
##                 r2 Pr(>r)
## Source      0.0000  1.000
## Sample_type 0.0257  0.859
## Vent_field  0.1259  0.818
## T_category  0.1359  0.547
## Substratum  0.0997  0.597
## Permutation: free
## Number of permutations: 999



ANOSIM tests on the clr-transformed data (Temp category and substratum).

## 
## Call:
## anosim(x = clr_euk9.dist, grouping = Substratum) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.2922 
##       Significance: 0.265 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.558 0.695 0.825 0.825 
## 
## Dissimilarity ranks between and within classes:
##         0%   25%  50%  75% 100%  N
## Between 11 17.25 21.5 26.5   33 14
## basalt  26 26.00 26.0 26.0   26  1
## sulfide  1  6.00 12.0 30.0   36 21
## 
## Call:
## anosim(x = clr_euk9.dist, grouping = T_category) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.05625 
##       Significance: 0.291 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.175 0.225 0.256 0.312 
## 
## Dissimilarity ranks between and within classes:
##         0%   25%  50%   75% 100%  N
## Between  4 12.75 18.5 24.25   36 20
## high     1  3.75  8.0 30.75   34 10
## low     11 17.00 24.5 27.50   29  6

RESULT: Neither temperature nor substratum affect variation in microeukaryote diveristy.


4.4.3 Archaea

Sample ECw11 produced insuffient archaeal reads so there are only 8 samples included here.


Envfit. No signficant predictors.

## 
## ***VECTORS
## 
##                   NMDS1       NMDS2     r2 Pr(>r)
## Temperature -7.3631e-05 -1.0000e+00 0.3679  0.303
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                             NMDS1    NMDS2
## SourceTubeworms            0.0000   0.0000
## Sample_typebk              4.1121   0.0005
## Sample_typed             -12.3363  -0.0015
## Vent_fieldClam_Bed       -12.3328  -0.0156
## Vent_fieldMain_Endeavour -12.3347   0.0012
## Vent_fieldMiddle_Valley   37.0032   0.0047
## T_categoryhigh           -12.3359  -0.0038
## T_categorylow             20.5598   0.0063
## Substratumbasalt         -12.3271   0.0094
## Substratumsulfide          1.7610  -0.0013
## 
## Goodness of fit:
##                 r2 Pr(>r)  
## Source      0.0000  1.000  
## Sample_type 0.0476  0.588  
## Vent_field  0.4286  0.266  
## T_category  0.2382  0.063 .
## Substratum  0.0204  0.890  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999


ANOSIM tests on the clr-transformed data (Temp category and substratum).

## 
## Call:
## anosim(x = clr_arch8.dist, grouping = Substratum) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: -0.483 
##       Significance: 1 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
##     1     1     1     1 
## 
## Dissimilarity ranks between and within classes:
##         0% 25% 50%  75% 100%  N
## Between  1   3   8 14.5   22  7
## basalt  NA  NA  NA   NA   NA  0
## sulfide  3  10  16 23.0   28 21
## 
## Call:
## anosim(x = clr_arch8.dist, grouping = T_category) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: -0.01538 
##       Significance: 0.476 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.282 0.364 0.405 0.569 
## 
## Dissimilarity ranks between and within classes:
##         0%   25%  50%   75% 100%  N
## Between  2  5.50 12.0 24.50   28 15
## high     7 10.75 14.5 18.25   21 10
## low      1 11.50 22.0 22.50   23  3

RESULT: Neither temperature or substratum have an effect on variation in archaeal diversity.


4.5 Microbes (all 23 samples)

4.5.1 Three domains, QPCR-balanced

Rel

Envfit. Temperature is the only significant predictor.

## 
## ***VECTORS
## 
##                NMDS1    NMDS2     r2 Pr(>r)
## Temperature  0.15956 -0.98719 0.0717  0.489
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                             NMDS1    NMDS2
## SourceBackground          -7.4027  29.4711
## SourceDiffuse              5.6560   3.7974
## SourceTubeworms           -1.3373  -8.8458
## Vent_fieldClam_Bed        26.0750  -8.7773
## Vent_fieldMain_Endeavour -34.6422  -1.8391
## Vent_fieldMiddle_Valley   73.4745  19.1430
## Substratumbasalt          20.8456  -3.7023
## Substratumnone            -7.4027  29.4711
## Substratumsulfide         -5.4680  -4.6601
## 
## Goodness of fit:
##                r2 Pr(>r)  
## Source     0.0282  0.902  
## Vent_field 0.2953  0.014 *
## Substratum 0.0395  0.819  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

ANOSIM tests (Vent field, substratum, sample type, temperature category)

## 
## Call:
## anosim(x = clr_balMic.dist, grouping = Vent_field) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.4485 
##       Significance: 0.002 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.169 0.225 0.263 0.319 
## 
## Dissimilarity ranks between and within classes:
##                 0%    25%   50%    75% 100%   N
## Between          1  99.75 155.0 209.75  253 154
## Clam_Bed         3  83.00 175.0 225.50  246  15
## Main_Endeavour   2  26.50  68.5 111.50  184  78
## Middle_Valley  120 177.50 180.0 187.75  204   6
## 
## Call:
## anosim(x = clr_balMic.dist, grouping = Source) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.2771 
##       Significance: 0.021 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.144 0.192 0.251 0.314 
## 
## Dissimilarity ranks between and within classes:
##            0%   25% 50%   75% 100%   N
## Between    13 93.50 136 186.0  252 151
## Background 36 43.00  50  58.0   66   3
## Diffuse     4 83.00 174 222.0  253  21
## Tubeworms   1 25.25  53 187.5  236  78
## 
## Call:
## anosim(x = clr_balMic.dist, grouping = Substratum) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.1286 
##       Significance: 0.155 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.178 0.242 0.302 0.389 
## 
## Dissimilarity ranks between and within classes:
##         0%    25%   50%    75% 100%   N
## Between  2  74.50 138.0 199.00  252 135
## basalt  15 189.75 230.5 246.25  253  10
## none    36  43.00  50.0  58.00   66   3
## sulfide  1  52.00 114.0 169.00  241 105
## 
## Call:
## anosim(x = clr_balMic.dist, grouping = T_category) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: -0.05103 
##       Significance: 0.584 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.185 0.276 0.347 0.407 
## 
## Dissimilarity ranks between and within classes:
##         0%   25%   50%    75% 100%   N
## Between  1 59.25 120.5 185.00  252 102
## high     6 18.50  26.0 207.50  236  15
## low      4 75.75 135.0 190.25  253 136

RESULT: Variation in microbial community diversity is not tied to any of the measured variables.


4.5.2 Bacteria

# Imputation of zeros.
codaBact23_cmultRepl <- cmultRepl(codaBact23nonsingleton, method = "CZM", 
    output = "p-counts")

# Need metadata for the next part. Must be in the same order as the
# samples.
metadata23 <- metadata23[match(rownames(codaBact23nonsingleton), rownames(metadata23)), 
    ]

# Centered log-ratio transformation
clr_codaBact23 <- apply(codaBact23_cmultRepl, 2, function(x) {
    log(x) - mean(log(x))
})

# Calculate Aitchison distance matrix
clr_codaBact23.dist <- dist(clr_codaBact23)

# NMDS
clr_codaBact23.NMDS <- metaMDS(clr_codaBact23.dist, k = 2)
clr_codaBact23.NMDS.df <- data.frame(x = clr_codaBact23.NMDS$points[, 1], 
    y = clr_codaBact23.NMDS$points[, 2], Source = as.factor(metadata23[, 
        1]), Vent_field = as.factor(metadata23[, 3]), Temperature = metadata23[, 
        4], T_category = as.factor(metadata23[, 5]), Substratum = as.factor(metadata23[, 
        6]))

ggplot(data = clr_codaBact23.NMDS.df, aes(x = x, y = y, color = Temperature, 
    shape = Source)) + geom_point(size = 4) + geom_text(aes(label = rownames(metadata23)), 
    hjust = 0.5, vjust = 1.5, size = 4) + ggtitle("Bacteria (23)") + scale_color_gradient(low = "blue", 
    high = "red")

Envfit. Temperature is the one and only significant predictor.

## 
## ***VECTORS
## 
##                NMDS1    NMDS2    r2 Pr(>r)    
## Temperature -0.83099  0.55628 0.485  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                             NMDS1    NMDS2
## Vent_fieldClam_Bed        -4.6834  -1.0941
## Vent_fieldMain_Endeavour   5.9497   7.3926
## Vent_fieldMiddle_Valley  -12.3115 -22.3849
## Substratumbasalt          37.8105  13.7544
## Substratumnone            -7.2824 -24.7856
## Substratumsulfide        -11.1470   0.3723
## SourceBackground          -7.2824 -24.7856
## SourceDiffuse             -6.7661 -20.9488
## SourceTubeworms            5.3238  16.9999
## 
## Goodness of fit:
##                r2 Pr(>r)
## Vent_field 0.0424  0.806
## Substratum 0.1294  0.221
## Source     0.1030  0.358
## Permutation: free
## Number of permutations: 999

ANOSIM tests (Vent field, substratum, sample type, temperature category)

## 
## Call:
## anosim(x = clr_codaBact23.dist, grouping = Vent_field) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: -0.1619 
##       Significance: 0.937 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.171 0.222 0.270 0.332 
## 
## Dissimilarity ranks between and within classes:
##                0%    25%   50%    75% 100%   N
## Between         3  58.25 117.5 177.75  247 154
## Clam_Bed        2  65.50  85.0 117.00  167  15
## Main_Endeavour  1 102.25 172.0 216.75  253  78
## Middle_Valley  31  43.00  60.5  85.50  109   6
## 
## Call:
## anosim(x = clr_codaBact23.dist, grouping = Substratum) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: -0.04545 
##       Significance: 0.563 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.206 0.268 0.324 0.423 
## 
## Dissimilarity ranks between and within classes:
##         0%   25% 50%    75% 100%   N
## Between  1 55.50 120 189.00  253 135
## basalt  32 86.25 166 195.75  241  10
## none     3  5.50   8  11.00   14   3
## sulfide  4 83.00 131 189.00  252 105
## 
## Call:
## anosim(x = clr_codaBact23.dist, grouping = Source) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: -0.07129 
##       Significance: 0.682 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.154 0.206 0.267 0.323 
## 
## Dissimilarity ranks between and within classes:
##            0%   25%   50%   75% 100%   N
## Between     2  68.5 115.0 170.5  247 151
## Background  3   5.5   8.0  11.0   14   3
## Diffuse     1   9.0  18.0  32.0   53  21
## Tubeworms  40 123.5 184.5 209.5  253  78
## 
## Call:
## anosim(x = clr_codaBact23.dist, grouping = T_category) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.4051 
##       Significance: 0.011 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.210 0.281 0.333 0.398 
## 
## Dissimilarity ranks between and within classes:
##         0%    25%   50%    75% 100%   N
## Between 70 110.75 148.5 207.25  253 102
## high    40  54.00 119.0 186.50  233  15
## low      1  34.75  79.0 174.25  247 136

RESULT: Variation in bacterial diversity is tied to temperature, but not to vent field, substratum or sample type.


4.5.3 Eukarya

Envfit. Sample type is a weak predictor.

## 
## ***VECTORS
## 
##               NMDS1   NMDS2     r2 Pr(>r)
## Temperature 0.12405 0.99228 0.0231  0.803
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                             NMDS1    NMDS2
## Vent_fieldClam_Bed         6.8767   4.6585
## Vent_fieldMain_Endeavour   8.3035  -5.6479
## Vent_fieldMiddle_Valley  -37.3013  11.3678
## Substratumbasalt          20.3966  -4.5481
## Substratumnone           -20.2410  -6.4232
## Substratumsulfide         -2.7507   2.8007
## SourceBackground         -20.2410  -6.4232
## SourceDiffuse            -25.4307 -32.3771
## SourceTubeworms           18.3645  18.9161
## 
## Goodness of fit:
##                r2 Pr(>r)  
## Vent_field 0.0872  0.446  
## Substratum 0.0421  0.771  
## Source     0.2486  0.022 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999


ANOSIM tests (Vent field, substratum, sample type, temperature category)

## 
## Call:
## anosim(x = clr_codaEuk23.dist, grouping = Vent_field) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: -0.04486 
##       Significance: 0.617 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.156 0.223 0.254 0.309 
## 
## Dissimilarity ranks between and within classes:
##                0%    25%   50%    75% 100%   N
## Between         1  62.25 123.0 192.75  253 154
## Clam_Bed        5  39.00  86.0 168.50  222  15
## Main_Endeavour  2  75.50 145.0 194.75  252  78
## Middle_Valley  95 100.50 108.5 136.00  199   6
## 
## Call:
## anosim(x = clr_codaEuk23.dist, grouping = Substratum) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.03239 
##       Significance: 0.379 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.197 0.247 0.314 0.374 
## 
## Dissimilarity ranks between and within classes:
##         0%   25% 50%   75% 100%   N
## Between  4  71.0 124 192.0  253 135
## basalt  14 154.5 192 227.5  252  10
## none     1  26.5  52  57.5   63   3
## sulfide  2  59.0 127 185.0  247 105
## 
## Call:
## anosim(x = clr_codaEuk23.dist, grouping = Source) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.1153 
##       Significance: 0.15 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.149 0.208 0.245 0.324 
## 
## Dissimilarity ranks between and within classes:
##            0%   25% 50%    75% 100%   N
## Between     4 77.50 131 192.50  252 151
## Background  1 26.50  52  57.50   63   3
## Diffuse     8 18.00  38  65.00  166  21
## Tubeworms   2 79.25 148 200.75  253  78
## 
## Call:
## anosim(x = clr_codaEuk23.dist, grouping = T_category) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: -0.1266 
##       Significance: 0.818 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.185 0.249 0.289 0.344 
## 
## Dissimilarity ranks between and within classes:
##         0%   25%   50%    75% 100%   N
## Between 12 67.25 102.5 159.50  247 102
## high     2  8.50  27.0 167.50  221  15
## low      1 81.75 147.5 198.25  253 136

RESULT: Variation in microeukaryote diversity is not tied to any of the measured variables.


4.5.4 Archaea


Envfit. Temperature and sample type are significant predictors.

## 
## ***VECTORS
## 
##               NMDS1   NMDS2    r2 Pr(>r)    
## Temperature 0.68773 0.72597 0.603  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                             NMDS1    NMDS2
## Vent_fieldClam_Bed         3.5115  -8.8294
## Vent_fieldMain_Endeavour  -4.3847   7.6991
## Vent_fieldMiddle_Valley    9.8609 -13.9853
## Substratumbasalt         -11.0887  -8.9927
## Substratumnone           -18.6048   0.7761
## Substratumsulfide          6.6780   2.2428
## SourceBackground         -18.6048   0.7761
## SourceDiffuse            -30.7007  -7.6776
## SourceTubeworms           22.5600   4.2846
## 
## Goodness of fit:
##                r2 Pr(>r)    
## Vent_field 0.0820  0.487    
## Substratum 0.0806  0.512    
## Source     0.4460  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999


ANOSIM tests (Vent field, substratum, sample type, temperature category)

## 
## Call:
## anosim(x = clr_codaArch22.dist, grouping = Vent_field) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.06507 
##       Significance: 0.297 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.170 0.213 0.258 0.312 
## 
## Dissimilarity ranks between and within classes:
##                0%   25%   50%    75% 100%   N
## Between         7 58.00 126.0 179.00  231 137
## Clam_Bed       22 53.75 156.5 197.25  228  10
## Main_Endeavour  1 60.75 103.5 156.75  226  78
## Middle_Valley  10 69.75 121.0 172.25  178   6
## 
## Call:
## anosim(x = clr_codaArch22.dist, grouping = Source) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.2482 
##       Significance: 0.02 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.157 0.200 0.232 0.279 
## 
## Dissimilarity ranks between and within classes:
##            0%   25%   50%    75% 100%   N
## Between     9 62.00 136.0 186.00  231 141
## Background  7  7.50   8.0  10.50   13   3
## Diffuse    14 52.00  79.0  99.00  167  21
## Tubeworms   1 63.75 113.5 152.75  218  66
## 
## Call:
## anosim(x = clr_codaArch22.dist, grouping = Substratum) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: -0.2803 
##       Significance: 0.985 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.198 0.269 0.326 0.360 
## 
## Dissimilarity ranks between and within classes:
##         0%  25% 50%   75% 100%   N
## Between  1 44.0  90 155.0  227 117
## basalt   4 66.0  81  93.0  111   6
## none     7  7.5   8  10.5   13   3
## sulfide  5 94.0 146 194.0  231 105
## 
## Call:
## anosim(x = clr_codaArch22.dist, grouping = T_category) 
## Dissimilarity: euclidean 
## 
## ANOSIM statistic R: 0.4778 
##       Significance: 0.001 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.167 0.228 0.262 0.340 
## 
## Dissimilarity ranks between and within classes:
##         0%    25%   50%    75% 100%   N
## Between 23 112.75 157.0 192.25  231  96
## high    45 103.50 134.0 156.50  218  15
## low      1  34.75  74.5 138.50  230 120

RESULT: Variation in archaeal diversity is tied to temperature and sample type.


5 Relative Abundance plots

Convert counts to relatvie abundances for stacked bar plots.

5.1 Macrofauna

Group “Others” include taxa that are always less than 2% relative abundance.

# Convert count data to relative abundances (it must be as a matrix)
macro9matrix <- as.matrix(macro9)
Macro9rel <- as.data.frame(make_relative(macro9matrix))

# Create row with maximum relative abundance for each species
colMax <- function(Macro9rel) sapply(Macro9rel, max, na.rm = TRUE)
MaxRel <- colMax(Macro9rel)
Macro9rel <- rbind(Macro9rel, MaxRel = MaxRel)

# Order columns based on the maximum relative abundance
Macro9rel <- Macro9rel[, order(-Macro9rel[nrow(Macro9rel), ])]

# Create 'Other' column with species always less than 2% relative
# abundance
Macro9rel$Others <- rowSums(Macro9rel[, c(17:41)])
MacTax <- Macro9rel[1:9, -c(17:41)]
MacTax$Sample <- rownames(MacTax)

# Create long table
MacTaxLong <- melt(MacTax, id = "Sample")

# Display taxa and samples in the desired order
MacTaxLong$variable <- factor(MacTaxLong$variable, levels = c("Lepetodrilus_fucensis", 
    "Depressigyra_globulus", "Paralvinella_palmiformis", "Amphisamytha_carldarei", 
    "Paralvinella_sulfincola", "Provanna_variabilis", "Hydrozoa_non_ID", 
    "Euphilomedes_climax", "Protomystides_verenae", "Fucaria_striata", 
    "Munidopsis_alvisca", "Branchinotogluma_tunnicliffeae", "Pardalisca_endeavouri", 
    "Helicoradomenia_juani", "Clypeosectus_curvus", "Paralvinella_pandorae", 
    "Others"))

MacTaxLong$Sample <- factor(MacTaxLong$Sample, levels = c("MVw13", "ECw11", 
    "EMw3", "EMw5", "MVw12", "ECw10", "EMw1", "EMw4", "EMw8"))

MACROcolors17 <- c("#fcd0a1", "#ad2929", "#0fa3b1", "#81f5f5", "#ffea51", 
    "#d68fd6", "#6144bf", "#f6ae2d", "#8c1c13", "#53917e", "#f26326", "#fffcc9", 
    "#7098d1", "#cfea9f", "#d26280", "#0686d6", "#979898")

ggplot(MacTaxLong, aes(x = Sample, y = value, fill = variable, order = Sample)) + 
    geom_bar(stat = "identity", width = 0.9) + scale_fill_manual(values = MACROcolors17, 
    label = c("L. fucensis", "D. globulus", "P. palmiformis", "A. carldarei", 
        "P. sulfincola", "P. variabilis", "Hydrozoa non ID", "E. climax", 
        "P. verenae", "F. striata", "M. alvisca", "B. tunnicliffeae", "P. endeavouri", 
        "H. juani", "C. curvus", "P. pandorae", "Others")) + theme(axis.text.x = element_text(angle = 90, 
    hjust = 1, vjust = 0.5), legend.text = element_text(size = 12), legend.title = element_blank()) + 
    scale_y_continuous(expand = c(0, 0)) + labs(x = NULL, y = "Relative Abundance")

5.2 Meiofauna

5.3 Bacteria

For all microbial plots, the group “Others” includes taxa that are always less than 1% relative abundance. Microbial OTUs were aggregated by high level taxonomies in Excel.

TaxaB23 <- read.csv("Bact23_high_level_taxa.csv", header = TRUE, sep = ",")
TaxaB23LONG <- melt(TaxaB23, id = c("Sample", "Type"))

TaxaB23LONG$variable <- factor(TaxaB23LONG$variable, levels = c("Alphaproteobacteria", 
    "Gammaproteobacteria", "Deltaproteobacteria", "Proteobacteria", "Bacteroidetes", 
    "Epsilonbacteraeota", "Acidobacteria", "Actinobacteria", "Aquificae", 
    "Chloroflexi", "Firmicutes", "Nitrospinae", "Patescibacteria", "Others", 
    "Bacteria"))

# Worms ordered low to high basal temperature
TaxaB23LONG$Sample <- factor(TaxaB23LONG$Sample, levels = c("bkEM", "bkEC", 
    "bkMV", "EMf3", "ECf11", "EMf2", "ECf9", "MVf12", "EMf7", "EMf1", "EMw6", 
    "EMw3", "ECw11", "EMw2", "EMw5", "ECw9", "MVw13", "MVw12", "ECw10", 
    "EMw7", "EMw1", "EMw4", "EMw8"))

ChromatBact15 <- c("#c60609", "#f3722c", "#f8961e", "#f9c74f", "#9fc681", 
    "#47682f", "#6bc4a9", "#327e67", "#91a9bd", "#577590", "#b1a0f1", "#7a5ee5", 
    "#5dbdd5", "#695746", "#928477")

labels <- c(Bkgd = "Bkgd", Diffuse = "Diffuse", Tubeworms = "Tubeworm grabs")

ggplot(TaxaB23LONG, aes(x = Sample, y = value, fill = variable, order = variable)) + 
    geom_bar(stat = "identity", width = 0.9) + facet_grid(. ~ Type, labeller = labeller(Type = labels), 
    scales = "free", space = "free") + scale_fill_manual(values = ChromatBact15, 
    name = "Bacteria", label = c("Alphaproteobacteria", "Gammaproteobacteria", 
        "Deltaproteobacteria", "Other Proteobacteria", "Bacteroidetes", 
        "Epsilonbacteraeota", "Acidobacteria", "Actinobacteria", "Aquificae", 
        "Chloroflexi", "Firmicutes", "Nitrospinae", "Patescibacteria", 
        "Others", "Unclassified Bacteria")) + theme(axis.text.x = element_text(angle = 90, 
    hjust = 1, vjust = 0.5), legend.text = element_text(size = 9), legend.title = element_text(size = 12)) + 
    scale_y_continuous(expand = c(0, 0)) + labs(x = NULL, y = "Relative Adundance")

5.4 Eukarya

TaxaE23 <- read.csv("Euk23_high_level_taxa.csv", header = TRUE, sep = ",")

5.5 Archaea

TaxaA22 <- read.csv("Arch22_high_level_taxa.csv", header = TRUE, sep = ",")

5.6 Balanced microbial community

MG <- read.csv("Combined_micro_major_groups.csv", header = TRUE, sep = ",")

MG$Taxa <- factor(MG$Taxa, levels = c("Bacteria-Gammaproteobacteria", "Bacteria-Alphaproteobacteria", 
    "Bacteria-Bacteroidetes", "Bacteria-Epsilonbacteraeota", "Bacteria-Deltaproteobacteria", 
    "Archaea-Thaumarchaeota", "Archaea-Euryarchaeota", "Bacteria-Patescibacteria", 
    "Bacteria-Chloroflexi", "Eukaryota-Alveolata", "Eukaryota-Rhizaria", 
    "Bacteria", "Bacteria-Actinobacteria", "Bacteria-Proteobacteria", "Bacteria-Betaproteobacteriales", 
    "Eukaryota-Opisthokonta", "Eukaryota", "Archaea-Asgardaeota", "Bacteria-Acidobacteria", 
    "Eukaryota-Amoebozoa", "Bacteria-Aquificae", "Archaea-Crenarchaeota", 
    "Bacteria-Other", "Archaea-Other", "Eukaryota-Other"))

# Ordered according to hierarchichal clustering dendrogram
MG$source <- factor(MG$source, levels = c("EMf11", "bkMV", "bkEC", "bkEM", 
    "EMf3", "EMf2", "ECf9", "EMf1", "EMf7", "MVf12", "MVw12", "EMw1", "ECw10", 
    "EMw7", "EMw4", "EMw8", "MVw13", "ECw11", "ECw9", "EMw3", "EMw6", "EMw2", 
    "EMw5"))

TAXAcolors25 <- c("#2bd9fe", "#fffcc9", "#45cb85", "#ad2929", "#46237a", 
    "#f26326", "#0686d6", "#fcd0a1", "#53917e", "#470ff4", "#cc2084", "#63535b", 
    "#e1dd8f", "#37a4bf", "#cfea9f", "#750d37", "#6144bf", "#7098d1", "#f6ae2d", 
    "#81f5f5", "#5c5fe8", "#f49390", "#979898", "#a0acad", "#bcb9b8")

ggplot(MG, aes(x = source, y = relabund, fill = Taxa, order = source)) + 
    geom_bar(stat = "identity", width = 0.9) + scale_fill_manual(values = TAXAcolors25, 
    labels = c("Gammaproteobacteria", "Alphaproteobacteria", "Bacteroidetes", 
        "Epsilonbacteraeota", "Deltaproteobacteria", "Thaumarchaeota", 
        "Euryarchaeota", "Patescibacteria", "Chloroflexi", "Alveolata", 
        "Rhizaria", "Uncl. Bacteria", "Actinobacteria", "Proteobacteria", 
        "Betaproteobacteria", "Opisthokonta", "Uncl. Eukaryota", "Asgardaeota", 
        "Acidobacteria", "Amoebozoa", "Aquificae", "Crenarchaeota", "Other Bacteria", 
        "Other Archaea", "Other Eukarya")) + theme(axis.text.x = element_text(angle = 90, 
    hjust = 1, vjust = 0.5), legend.title = element_text(size = 14, face = "bold"), 
    legend.text = element_text(size = 9)) + scale_y_continuous(expand = c(0, 
    0)) + labs(x = NULL, y = "Relative Abundance")

5.7 Details of top 4 bacterial groups

Code only shown for Alphaproteobacteria.

alphas <- read.csv("Alphas.csv", header = TRUE, sep = ",")
alphasLONG <- melt(alphas, id = c("Sample"))

AlphaColors <- c("#cc2084", "#63535b", "#e1dd8f", "#37a4bf", "#cfea9f", 
    "#750d37", "#6144bf", "#7098d1", "#f6ae2d", "#81f5f5", "#5c5fe8", "#f49390", 
    "#934b00", "#ffea51", "#488243", "#fffcc9", "#45cb85", "#ad2929", "#46237a", 
    "#f26326", "#0686d6", "#fcd0a1", "#53917e", "#470ff4", "#CCC9CF", "#CCC9CF", 
    "#CCC9CF")

alphasLONG$variable <- factor(alphasLONG$variable, levels = c("Rhodobacteraceae", 
    "SAR11_clade", "unknown_Alphaproteobacteria", "Roseobacter_clade_NAC11_7_lineage", 
    "Parvibaculales_PS1_clade", "Kordiimonadales", "Robiginitomaculum", 
    "Euryhalocaulis", "Sedimentitalea", "Bradyrhizobium", "Pseudahrensia", 
    "Rhodospirillales_AEGEAN_169_marine_group", "Devosiaceae", "Rhizobiales", 
    "Hellea", "Methylobacterium", "Beijerinckiaceae", "Aestuariispira", 
    "Nisaeaceae_OM75_clade", "Kiloniellaceae", "Parvularculaceae", "Magnetospiraceae", 
    "Parvibaculaceae", "Other_Alphaproteobacteria"))

alphasLONG$Sample <- factor(alphasLONG$Sample, levels = c("ECf11", "bkMV", 
    "bkEC", "bkEM", "EMf3", "EMf2", "ECf9", "EMf1", "EMf7", "MVf12", "MVw12", 
    "EMw1", "ECw10", "EMw7", "EMw4", "EMw8", "MVw13", "ECw11", "ECw9", 
    "EMw3", "EMw6", "EMw2", "EMw5"))

ggplot(alphasLONG, aes(x = Sample, y = value, fill = variable, order = variable)) + 
    geom_bar(stat = "identity", width = 0.9) + scale_fill_manual(values = AlphaColors, 
    name = NULL, labels = c("Rhodobacteraceae", "SAR11 clade", "Unknown Alphaproteobacteria", 
        "Roseobacter clade NAC11-7", "Parvibaculales clade PS1", "Kordiimonadales", 
        "Robiginitomaculum", "Euryhalocaulis", "Sedimentitalea", "Bradyrhizobium", 
        "Pseudahrensia", "Rhodospirillales clade AEGEAN-169", "Devosiaceae", 
        "Rhizobiales", "Hellea", "Methylobacterium", "Beijerinckiaceae", 
        "Aestuariispira", "Nisaeaceae clade OM75", "Kiloniellaceae", "Parvularculaceae", 
        "Magnetospiraceae", "Parvibaculaceae", "Other Alphaproteobacteria")) + 
    theme(axis.text.x = element_text(size = 11, angle = 90, hjust = 1, 
        vjust = 0.5), legend.text = element_text(size = 9)) + labs(x = NULL, 
    y = "Relative Abundance")

gammas <- read.csv("Gammas.csv", header = TRUE, sep = ",")

bacteroidetes <- read.csv("Bacteroidetes.csv", header = TRUE, sep = ",")

epsilons <- read.csv("Epsilons.csv", header = TRUE, sep = ",")


6 Enriched species/OTUs

Species and OTUs that differentiated between different sample types (tubeworm grabs, diffuse fluids and background fluids) and grab sample temperature categories (highT, lowT) were identified using ALDEx2. Microbes were analyzed from all 13 grab samples (6 highT, 7 lowT) and 10 fluid samples (7 diffuse, 3 background) using both individual domains and the QPCR-balanced microbial assemblage. Fauna were analyzed from 9 grab samples (5 highT, 4 lowT). We include an analysis example below, but do not include all tests for brevity.

library(ALDEx2)

# create the 2-state condition for comparison
tempx2 <- as.vector(metadata9$T_category)
tempx2
## [1] "low"  "low"  "high" "high" "high" "high" "low"  "high" "low"
# ALDEx2 modular mode
flip_mac9 <- as.data.frame(t(macro9))
clr_mac9byT <- aldex.clr(flip_mac9, tempx2, mc.samples = 128, denom = "all")
kw_mac9byT <- aldex.kw(clr_mac9byT)
kw_mac9byT.effect <- aldex.effect(clr_mac9byT, tempx2, verbose = FALSE)
mac9byT.all <- data.frame(kw_mac9byT, kw_mac9byT.effect)

# mac9byT.all has details of differential taxa between highT and lowT
# (plotted in red below)

aldex.plot(mac9byT.all, type = "MW", test = "glm", xlab = "Dispersion", 
    ylab = "Difference", all.cex = 1.5, all.pch = 16, all.col = "lightgrey", 
    rare.col = "black")

  • In this case, macrofauna enriched in highT assemblages (effect < -1, Benjamini-Hochberg corrected p-values < 0.1 for GLM test):
    • Paralvinella sulfincola
    • Paralvinella pandorae
    • Depressigyra globulus
    • Paralvinella palmiformis
    • Branchinotogluma tunnicliffeae
  • Macrofauna enriched in lowT assemblages (effect > 1):
    • Amphisamytha caldarei
    • Helicoradomenia juani


7 Plotting enriched species/OTUs

Microbes presented here are from ALDEx2 tests using the QPCR-balanced microbial assemblage. Taxa with absolute effect values >1 were compiled in Excel for plotting.

7.1 HighT vs. lowT grabs

htlteBAL <- read.csv("enriched_HTvsLT_grabs.csv", header = TRUE, sep = ",")
htlteBAL$range <- ifelse(htlteBAL$Effect < 0, "below", "above")  # 'above'= low T   'below'= high T
htlteBAL <- htlteBAL[order(htlteBAL$Effect), ]
htlteBAL$Identity <- factor(htlteBAL$Identity, levels = htlteBAL$Identity)  # keeps the order for the plot rather than going alphabetical

ggplot(htlteBAL, aes(x = Identity, y = Effect, label = Effect)) + geom_bar(stat = "identity", 
    aes(fill = Domain), width = 0.6) + scale_fill_manual(name = NULL, labels = c("Archaea", 
    "Bacteria", "Eukarya", "Macrofauna", "Meiofauna"), values = c(Archaea = "#abdafc", 
    Bacteria = "#0686d6", Eukarya = "#cfea9f", Macrofauna = "#f6ae2d", 
    Meiofauna = "#f26326")) + theme(axis.text.x = element_text(size = 8), 
    axis.title.x = element_text(size = 9), axis.text.y = element_text(size = 6), 
    axis.title.y = element_blank(), legend.position = "none") + ggtitle("HighT vs. LowT grabs") + 
    scale_x_discrete(position = "bottom") + coord_flip()

7.2 HighT grabs vs. diffuse

All 6 highT grab samples and 3 associated diffuse fluid samples.

7.3 LowT grabs vs. diffuse

All 7 lowT grab samples and 4 associated diffuse fluid samples.


8 Congruence

Symmetric procrustes analysis (function ‘procrustes’) and significance testing (function ‘protest’) performed on pairs of NMDS ordinations.

8.1 Eight sample NMDS

Archaeal reads were insufficient from sample ECw11 so NMDS plots without this sample for the other size classes and microbial domains will be necessary for comparisons.

8.2 Nine sample, QPCR-balanced microbial NMDS

Envfit.

## 
## ***VECTORS
## 
##               NMDS1   NMDS2     r2 Pr(>r)  
## Temperature 0.91082 0.41280 0.4431  0.026 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##                             NMDS1    NMDS2
## SourceTubeworms            0.0000   0.0000
## Sample_typebk            -25.9926  -0.0094
## Sample_typed              51.9852   0.0189
## Vent_fieldClam_Bed        51.9647   0.0370
## Vent_fieldMain_Endeavour -41.5885   0.0105
## Vent_fieldMiddle_Valley   52.0066  -0.0633
## T_categoryhigh            51.9988  -0.0145
## T_categorylow            -64.9985   0.0181
## Substratumbasalt         -65.0383 -67.4684
## Substratumsulfide         18.5824  19.2767
## 
## Goodness of fit:
##                 r2 Pr(>r)  
## Source      0.0000  1.000  
## Sample_type 0.1000  0.469  
## Vent_field  0.1600  0.845  
## T_category  0.2502  0.025 *
## Substratum  0.1857  0.190  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999


8.3 Macrofauna vs. Meiofauna

library(vegan)

clr_mac9meio9NMDS <- procrustes(clr_mac9.NMDS, clr_meio9.NMDS, scale = FALSE, 
    choices = c(1, 2), symmetric = TRUE)

par(mfrow = c(1, 2))
plot(clr_mac9meio9NMDS, kind = 1)
plot(clr_mac9meio9NMDS, kind = 2)

Testing for significance.

protest(clr_mac9.NMDS, clr_meio9.NMDS, scores = "sites", permutations = how(nperm = 999))
## 
## Call:
## protest(X = clr_mac9.NMDS, Y = clr_meio9.NMDS, scores = "sites",      permutations = how(nperm = 999)) 
## 
## Procrustes Sum of Squares (m12 squared):        0.2588 
## Correlation in a symmetric Procrustes rotation: 0.8609 
## Significance:  0.005 
## 
## Permutation: free
## Number of permutations: 999

8.4 Macrofauna vs. Bacteria

## 
## Call:
## protest(X = clr_mac9.NMDS, Y = clr_bact9.NMDS, scores = "sites",      permutations = how(nperm = 999)) 
## 
## Procrustes Sum of Squares (m12 squared):        0.6513 
## Correlation in a symmetric Procrustes rotation: 0.5905 
## Significance:  0.096 
## 
## Permutation: free
## Number of permutations: 999

8.5 Macrofauna vs. Archaea

## 
## Call:
## protest(X = clr_mac8.NMDS, Y = clr_arch8.NMDS, scores = "sites",      permutations = how(nperm = 999)) 
## 
## Procrustes Sum of Squares (m12 squared):        0.6865 
## Correlation in a symmetric Procrustes rotation: 0.5599 
## Significance:  0.02 
## 
## Permutation: free
## Number of permutations: 999

8.6 Macrofauna vs. Eukarya

## 
## Call:
## protest(X = clr_mac9.NMDS, Y = clr_euk9.NMDS, scores = "sites",      permutations = how(nperm = 999)) 
## 
## Procrustes Sum of Squares (m12 squared):        0.9062 
## Correlation in a symmetric Procrustes rotation: 0.3062 
## Significance:  0.818 
## 
## Permutation: free
## Number of permutations: 999

8.7 Meiofauna vs. Bacteria

## 
## Call:
## protest(X = clr_meio9.NMDS, Y = clr_bact9.NMDS, scores = "sites",      permutations = how(nperm = 999)) 
## 
## Procrustes Sum of Squares (m12 squared):        0.6047 
## Correlation in a symmetric Procrustes rotation: 0.6287 
## Significance:  0.057 
## 
## Permutation: free
## Number of permutations: 999

8.8 Meiofauna vs. Archaea

## 
## Call:
## protest(X = clr_meio8.NMDS, Y = clr_arch8.NMDS, scores = "sites",      permutations = how(nperm = 999)) 
## 
## Procrustes Sum of Squares (m12 squared):        0.8721 
## Correlation in a symmetric Procrustes rotation: 0.3576 
## Significance:  0.285 
## 
## Permutation: free
## Number of permutations: 999

8.9 Meiofauma vs. Eukarya

## 
## Call:
## protest(X = clr_meio9.NMDS, Y = clr_euk9.NMDS, scores = "sites",      permutations = how(nperm = 999)) 
## 
## Procrustes Sum of Squares (m12 squared):        0.776 
## Correlation in a symmetric Procrustes rotation: 0.4733 
## Significance:  0.268 
## 
## Permutation: free
## Number of permutations: 999

8.10 Bacteria vs. Eukarya

## 
## Call:
## protest(X = clr_bact9.NMDS, Y = clr_euk9.NMDS, scores = "sites",      permutations = how(nperm = 999)) 
## 
## Procrustes Sum of Squares (m12 squared):        0.3398 
## Correlation in a symmetric Procrustes rotation: 0.8125 
## Significance:  0.001 
## 
## Permutation: free
## Number of permutations: 999

8.11 Archaea vs. Eukarya

## 
## Call:
## protest(X = clr_arch8.NMDS, Y = clr_euk8.NMDS, scores = "sites",      permutations = how(nperm = 999)) 
## 
## Procrustes Sum of Squares (m12 squared):        0.9796 
## Correlation in a symmetric Procrustes rotation: 0.1429 
## Significance:  0.378 
## 
## Permutation: free
## Number of permutations: 999

8.12 Bacteria vs. Archaea

## 
## Call:
## protest(X = clr_bact8.NMDS, Y = clr_arch8.NMDS, scores = "sites",      permutations = how(nperm = 999)) 
## 
## Procrustes Sum of Squares (m12 squared):        0.9622 
## Correlation in a symmetric Procrustes rotation: 0.1944 
## Significance:  0.756 
## 
## Permutation: free
## Number of permutations: 999

8.13 Microbes vs. Macrofauna

## 
## Call:
## protest(X = clr_mac9.NMDS, Y = clr_balMic9.NMDS, scores = "sites",      permutations = how(nperm = 999)) 
## 
## Procrustes Sum of Squares (m12 squared):        0.8275 
## Correlation in a symmetric Procrustes rotation: 0.4153 
## Significance:  0.445 
## 
## Permutation: free
## Number of permutations: 999

8.14 Microbes vs. Meiofauna

## 
## Call:
## protest(X = clr_meio9.NMDS, Y = clr_balMic9.NMDS, scores = "sites",      permutations = how(nperm = 999)) 
## 
## Procrustes Sum of Squares (m12 squared):        0.7268 
## Correlation in a symmetric Procrustes rotation: 0.5226 
## Significance:  0.093 
## 
## Permutation: free
## Number of permutations: 999


9 Proportionality

Microbial OTUs were binned by shared taxonomic assignment into 719 unique taxonomic identities. Prior to proportionality, rare taxa (i.e. those with counts of <2 in 7+ samples) were removed to minimize false discovery rate, resulting in final numbers of taxa for analysis of 335 unique microbial taxa, 14 macrofauna, and 17 meiofauna. Only rho values >0.75 are used in Cytoscape network mapping.

library(propr)
NoRares <- read.csv("Non_rares_for_propr.csv", header = T, row.names = 1)

# Imputation of zeros
NoRares_cmultRepl <- cmultRepl(NoRares, method = "CZM", output = "p-counts")
## No. corrected values:  85
# Centered log-ratio transformation
clr_NoRares <- apply(NoRares_cmultRepl, 2, function(x) {
    log(x) - mean(log(x))
})

# Proportionality analysis
NoRares.propr <- propr(counts = NoRares_cmultRepl, metric = "rho", p = 1000)
NoRaresCorrALL <- getResults(NoRares.propr)

10 Session Info

sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] ALDEx2_1.20.0       propr_4.2.6         zCompositions_1.3.4
##  [4] truncnorm_1.0-8     NADA_1.6-1.1        survival_3.1-12    
##  [7] MASS_7.3-51.6       compositions_2.0-0  funrar_1.4.1       
## [10] ggpubr_0.4.0        ggrepel_0.8.2       ggplot2_3.3.2      
## [13] reshape2_1.4.4      vegan_2.5-6         lattice_0.20-41    
## [16] permute_0.9-5       knitr_1.29         
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-148                bitops_1.0-6               
##  [3] matrixStats_0.56.0          GenomeInfoDb_1.24.2        
##  [5] tensorA_0.36.1              tools_4.0.2                
##  [7] backports_1.1.9             R6_2.4.1                   
##  [9] BiocGenerics_0.34.0         mgcv_1.8-31                
## [11] colorspace_1.4-1            withr_2.2.0                
## [13] gridExtra_2.3               tidyselect_1.1.0           
## [15] bayesm_3.1-4                curl_4.3                   
## [17] compiler_4.0.2              Biobase_2.48.0             
## [19] formatR_1.7                 DelayedArray_0.14.1        
## [21] labeling_0.3                scales_1.1.1               
## [23] DEoptimR_1.0-8              robustbase_0.93-6          
## [25] stringr_1.4.0               digest_0.6.25              
## [27] foreign_0.8-80              rmarkdown_2.3              
## [29] rio_0.5.16                  XVector_0.28.0             
## [31] pkgconfig_2.0.3             htmltools_0.5.0            
## [33] rlang_0.4.7                 readxl_1.3.1               
## [35] generics_0.0.2              farver_2.0.3               
## [37] BiocParallel_1.22.0         dplyr_1.0.2                
## [39] zip_2.1.1                   car_3.0-9                  
## [41] RCurl_1.98-1.2              magrittr_1.5               
## [43] GenomeInfoDbData_1.2.3      Matrix_1.2-18              
## [45] Rcpp_1.0.5                  munsell_0.5.0              
## [47] S4Vectors_0.26.1            abind_1.4-5                
## [49] lifecycle_0.2.0             stringi_1.5.3              
## [51] yaml_2.2.1                  carData_3.0-4              
## [53] SummarizedExperiment_1.18.2 zlibbioc_1.34.0            
## [55] plyr_1.8.6                  grid_4.0.2                 
## [57] parallel_4.0.2              forcats_0.5.0              
## [59] crayon_1.3.4                haven_2.3.1                
## [61] cowplot_1.1.0               splines_4.0.2              
## [63] hms_0.5.3                   pillar_1.4.6               
## [65] GenomicRanges_1.40.0        ggsignif_0.6.0             
## [67] stats4_4.0.2                glue_1.4.2                 
## [69] evaluate_0.14               data.table_1.13.0          
## [71] vctrs_0.3.4                 cellranger_1.1.0           
## [73] gtable_0.3.0                purrr_0.3.4                
## [75] tidyr_1.1.2                 xfun_0.17                  
## [77] openxlsx_4.1.5              broom_0.7.0                
## [79] rstatix_0.6.0               tibble_3.0.3               
## [81] IRanges_2.22.2              cluster_2.1.0              
## [83] ellipsis_0.3.1